Optics in the Schwarzschild space-time 
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Abstract 

Realistic modelling of radiation transfer in and from variable accretion disks around black 
holes requires the solution of the problem: find the constants of motion and equation of 
motion of a light-like geodesic connecting two arbitrary points in space. Here we give the 
complete solution of this problem in the Schwarzschild space-time. 



1 Introduction 

The first light detected from a relativistic region about a black hole was discovered by the ASCA 
satellite ( Tanaka et al. 199,4 Yaqoob et al. 1998h . The now accepted theoretical model describing 

the broad X-ray emission lines is that of an accretion disk around either a Kerr or a Schwarzschild 

black hole (iRevnolds et"aT1ll99,4lTurner et alll"997ilNandra et al.lll997ilFanton et al.lll997LlBromlev et al 
M£i iDabrowski Lasenbvll200ll ICadez et al.lll998l iMartocchia et al.ll2000h. This discovery i 



in- 



creased the interest for phenomena occurring in the vicinity of black holes. We now know that other 



intere sting high energy phenomena, such as X-ray flares feagano: 



2003) or quasi-periodic oscillations ( Strohmaver 2001 , Miller et al. 2002^ 1 are occuring in the en- 



et al .M200 1L iGoldwurm et al 



vironment of black holes. To further the understanding of such phenomena from the theoretical 
point of view, it is important to develop tools to model the phenomena themselves as well as to 
model radiation transfer in and from these strongly curved regions of space-time. 

Most current accretion disk models are very simple from the radiation transfer point of view. 
They consider disks as geometrically thin and optically thick. Light that reaches a far observer, 
comes directly (without beeing scattered) from a very definite point on the disk. Therefore, line 
profiles can be calculated by aiming light-like goedesic from a point on the disk to the observer 
or vice versa, aiming from the observer to points on the disk. In suc h ray-tracing procedures 
geodesic equations are usually so lved by direct numerical integration (jLaorJh99ll iPariev et al 



2001L IRevnolds &: Begelmanlll997f ). However, when modelling transient phenomena produced by 
small debris around black holes, or by other transient phenomena in the disk (moving hot spots, 
varying external illumination, waves), it becomes necessary to solve a more difficult radiation 
transfer problem, the problem of following a single photon through its more t han one scatterin 
and/or more than one possible path from the source to the eye of the observer (|Schnittmanll200 
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Schnittman fc Bertschingerl 12004 ) . The problem is further complicated by noting that photon 



arrival times from the same initial source to the observer may (and often will) be markedly different 
for photons reaching the observer along different possible paths. It is obvious, that a multiple 
scattering path cannot be effectively constructed by aiming geodesies between succesive scattering 
points as the number of succesive iterations required would soon blow up. The solution to the 
radiation transfer problem thus requires one to be able to follow a light ray from one to the other 
scattering point along its path. Thus, one would like to find the quickest way to determine all 
constants of motion of a geodesic that connects the given initial and final point. In this article 
we give analytic expressions which completely solve this problem in the Schwarzschild space-time. 
We also turn this analytic tool into a numerical code and demonstrate that it is much faster, 
more accurate and more transparent than aiming and integrating geodesic equations. This tool 
thus opens the possibility to solve complex radiation transfer problems in curved space-time using 
similar Monte Carlo techniques that are used in solving r a diatio n pro blems in flat space-t i me. 

Our work starts with the results of Chandrasekhar l|l992ft and iRauch fc Blandfordl (|l994h . 
who expressed the solutions to geodesic e quations in terms of elliptic integrals. By inverting their 



expressions into Jacobi elliptic functions ( Cadez et al. 1998h . we obtain simple solutions for all the 



three types of orbit equations that occur in Schwarzscild space-time. These solutions no longer 
contain branch ambiguities. Since these orbits are essentialy planar, their equations depend only on 
two nontrivial constants of motion: the angular momentum and the longitude of the periastron. 
Expressing the orbit equation at the initial and final point on the geodesic, one obtains two 
nonlinear equations for the two nontrivial constants of motion. However, since the longitude of 
the periastron occurs only linearly as the argument of elliptic functions, it is possible to use the 
elliptic functions addition theorem to elliminate the longitude of the periastron and obtain a single 
nonlinear equation for the angular momentum as function of initial and final coordinates. Here we 
derive these equations for all the three types of orbits and discuss their properties and solutions. 
We also write down all the other constants of motion in terms of final and initial point coordinates 
and give analytic expressions for travel times. 



2 Connecting two points with a light-like geodesic 



In the Schwarzschild space-time it is customary to introduce Schwarzscild coordinates t, r, 9 and 
if. In these coordinates geodesies are governed by the Hamiltonian: 



H 



1 



1 



2A1 



1 - 



2M 



1 

+ ~K 



1 



sin 



e 



IK 



which admits 8 constants of motion: the value of the Hamiltonian (H) and the value of the 
Lagrangian (L) (the ratio of the two defining the relation between time and proper time), the 
energy E = p t , the angular momentum (/), the longitude of the periastron (u) and the time of 
periastron passage. 

The angular momentum is expressed as / = / • n, where h is the unit vector along /, which is 
defined by the inclination of the orbit (e) and the longitude of the ascending node (Q). 

Consider a light-like geodesic connecting the initial point 3^ and the final point CP/. Five 
constants H — 0, L — 0, E, Q and e are determined readily. So one can use the true anomaly A 
as a parameter along the geodesic. From the initial to the final point A increases by: 



A A/ = 5X + 27rk , 



(2) 
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Figure 1: Left: The orbital plane in equatorial coordinates: n unit normal, e inclination, f2 longitude of 
the ascending node, uj longitude of the periastron and A the true anomaly. Right: The initial !Pj = (ft) 
and the final 7 f = (Of,(ff) point. 

where 

5\ = arccos [cos9iCos9 f + sin9ism 9 f cos(iff — (pi)~\ , (3) 

and the winding number k = ... — 1,0,1... tells how many times a geodesic winds around the black 
hole. 1 

The two constants Q and e are obtained from angular coordinates of the initial and final point 
(see Figure with the help of basic spherical trigonometry: 

sin 9; sin 9f sm({p f — ipA 

cose = { — ^ 4 

smoA 

Q sin 9i cos 9 1 sin oja — cos 9 t sin 9 r sin uj t . . 

-^g^-Q — — i i _± (O ) 

2 sine sin 5\ + sin 0, cos cos <£>j — cos ^ sin #j cosy?/ 
Since the variable A is conjugate to the orbital angular momentum, it obeys t he Poisson bracket 



relat ion [A, I] = 1, so that the angles 9 and <p are expressed with A as follows (ICadez &: Gomboc 
19961 eqs. (6 - 16)): 



cos 6* = — sinesin(A + uj) (6) 

(p — Q cose sin(A + uj) . . 

tan — — — it) 

2 sin 9 + cos(A + uj) 



and the differential equation for orbits of light-like geodesies becomes: 

du 



— = ±^a?-v?{l-u) . (8) 

Here u = 2M/r and a = 2ME/1. The solutions to this equation, called orbit equations, depend 
only on the parameter a and are of three types (Figure |2J): 

- type A: scattering orbits with both endpoints at infinity; their angular momentum param- 
eter is on the interval < a < Scattering orbits can never extend below r = 3M. 



^ote, if AA/ < 0, replace AA/ — > -AA/ and 



n — > —n. 



2 Note, that sin# = — cos 2 i 
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- type B: plunging orbits with one end at infinity and the other behind the horizon; a > 

- type C: near orbits with both ends behind the horizon of the black hole; their angular 
momentum parameter is on the interval < a < Near orbits can never reach beyond 
r = 3M. 

For completeness, we list the solutions of the orbit equation (JSJ) in terms of a few auxiliary param- 
eters. 

Type A: Introduce a = sin ^ (0 < ip < tt) and the following auxiliary parameters, func- 
tions of ip only: 

u x = i(l + 2 cos |) (9a) 

1 / i> — 27T\ . , , 

u 2 = - ( 1 + 2 cos — - — ) (9b) 



3 V 3 

M3 = i(l + 2co S i±^) (9c) 
m = (9d) 

Ml - U 3 

2 , s 

n = ; 9e 
Xi such that iti = u 2 — (u 2 — 1*3) cos 2 x« , (9f ) 

where U{ = 2M/rj and similarly Uf = 2M/rj. In terms of these, the type A differential equation 
for the orbit becomes (Figure |2l left): 

u = u 2 — (u 2 — ■u 3 )cn 2 (F(x i |m) + — Ira) , (10) 

v n ' 



with AA = A— Aj and AX/n e {x min — X{ , (2K(m)— x min )— x^}, where x min = F(arccos y/u 2 / (u 2 — tt 3 ) 
and Xj = F(arccos a/ (u 2 — ttj)/ (w 2 — u^)\m), se e Figure El left. Here F and K are the incomplete 
and complete elliptic integrals of the first kind (jWolfrarrl 1996J ) . 



Critical A: A special limiting case of solution (|TU|) for ijj = n (a = and U{ < 2/3 (r* > 3M) 
is: 

// = h tanh 2 (ArtanhA/u l + 1/3 H ) , (11) 

3 2 



with A A = A - Xi and A A G ^ArtanliA/I/^ - 2Artanh A /« i + 1/3 ,00}. 

Type B: Introduce a = ^^cosh/i (0 < /x < 00) and the auxiliary parameters, functions of 
fjL only: 

ui = -fl - 2 cosh ^) (12a) 



3 V 3 

m =±(l-—^Z}=\ (12b) 
2 y 2V«i(3«!-2) 1 

n= («i(3wi-2))"' (12c) 



1 1 — cos x« 

X« such that itj =U\-\ — ~ . (12d) 

n 2 1 + cos Xi 
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The type B differenial equation for the orbit is (Figure El centre): 



1 - cn(F(xi 


m) + ^ 


rnj 


l + cn(F{ Xi 


m) + ^ 


mj 



u = u i + -2 , , -;J L; . a Ai ...( . (13) 



with AA = A - Ai and AA/n G {F(xbh|^) - , F(xooM - F(xi\m)}, where ^sif = 

1— n 2 (l— ui) j l+n 2 ui 

arccos 1 , 2 ), — 4 and Xoo = arccos , , 1 . 

Type C : Introduce a = ^= sin ^ (0 < ^ < 7r) and the auxiliary parameters, functions of ip 
only: 

ui = i(l + 2cos|) (14a) 
u 2 = - l + 2cos— - — (14b) 



3 V 3 

U3 = 31 1 + 2C ° S 3 — ) (14C) 

m= 1 -(l- ^-^ + ^)/2 \ (Md) 

2 V V ( U l _ M 3)(«l - U 2 ) J 

n = ((u 1 -u 2 )(ui-u 3 ))~* (14e) 



Xi such that Uj = u\ H — — , (14f ) 

rr 1 + cos Xi 



The type C differential equation for the orbit is (Figure El right): 

1 1 -cn(F(YJm) + — \m) , 



with AA = A — Ai and AX/n G {— F(x_Bi?|m) — F(%j|m) , -F(xs#|m) — F(xi|m)}, where 



l-n 2 (l-ui) 

Xbh = arccos ^j^j . 



ft a 



3^ 



Critical C: A special limiting case of solution (|13|) for /i = (or (|15j) for ip 

and ^ > 2/3 (r< < 3M) is: 

w = + coth 2 (Arcothv/Mi + l/3 + ^) , (16) 

with A A = A - Ai and A A G {2Arcoth v /4/3 - 2ArcothA/Mi + 1/3 ,00}. 

The two constants of motion u and a can now be determined in two steps. First we determine 
the type of the orbit following reasoning illustrated in Figure El and t hen we solve the two equations 
obtained from the orbit equati on at the initial and the final point ( Brainik 19991 Gomboc 2001 
Cadez et al.ll2003L lKostidl2003h . 



We note that the parameter u appears in a contrived form as the starting point of the true 
anomaly only in the argument of the Jacobi functions, therefore, it can be eliminated from the 
two equations by using the Jacobi elliptic functions addition theorem. Let v = F(xi\m) + AA/n 
be the argument of elliptic functions at the final point of the orbit and z = F(xi\ m ) be the argu- 
ment of elliptic functions at the initial point of the orbit. Then one can use orbit equations to write: 
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Figure 2: Left: Some orbits of type A with parameters ri = 20M and a cr a — a G 

{0.2,0.15,0.1,0.05,0.005}. Middle: Orbits of type B with parameters r* = 20M and a - a crit G 

{2,0.7,0.2,0.05,0.005}. Right: Orbits of type C with parameters u = 2.00001M and a crit — a £ 
{0.01,0.005,0.001}. The radius of the black circle is the Schwarzschild radius. 




Figure 3: For a known initial point CP^ , the two critical orbits through this point divide the orbital plane 
into three regions. Left: initial point at r > 3M. If the final point CP/ is in region 3, then all orbits from 
CPj are of type A. If CP/ is in region 1, then all orbits from CPj are of type B. If CP y is in region 2, then orbits 
leading to Tf are of type A if the trajectory winds around the black hole, and of type B if it goes directly 
from 3^ to CP/. Right: initial point at r < 3M. If the final point 7 f is in region 3, then all orbits from 
are of type B. If Tj is in region 1, then orbits from "Pi are of type C. If 7 f is in region 2, then orbits 
leading to CP / are of type C if the trajectory winds around the black hole, and of type B if it goes directly 
from CPj to CP/. 
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Figure 4: Functions cn(x\m) and sn(x|m) along orbits: left type A, middle type B, right type C. x m i n 
and x max are at the endpoints of the orbit and are designated by 7- and if they are at spatial infinity, 
and Th- and < Pf l + if they are at the horizon of the black hole. Interval definitions for the endpoints follow 
after equations (|TT))). (|T31) and (fT31) . Xi = z and Xf = v are at the initial and the final point of the orbit. 



Type A: 



cn (v\m 



sn l vim 



dn 2 



v\m 



cn (z\m 



sn [z \m 



dn 2 ( 



z\m 



u 2 


-Uf 


u 2 


- u 3 


u f 


- u 3 


u 2 


- u 3 


Ul 


~ u f 


Ml 


- u 3 


u 2 


- Ui 


u 2 


- u 3 


Ui 


- u 3 


u 2 


- u 3 




- Ui 




- u 3 



Types B and C: 

1 — n 2 {uf — ui) 



(17a) 
(17b) 
(17c) 
(17d) 
(17e) 
(17f) 



cn (v 



sn 2 (v 



dn 2 



cn z 



sn 2 (2; 



dn 2 (^ 



m) 



m) 



m) 



m) 



m) 



m) 



1 + n 2 {uf — Mi) 
i.n 2 {uf — Ui) 
{l + n 2 {u f -ui)y 
Amn 2 {uf — Mi) 



(1 + n 2 (uf — Mi)) 2 
n 2 (ui - ui) 



1 + n 2 (ui — Ui) 
An 2 {ui — Mi) 
(1 +n 2 (ui -Mi)) ; 



1 - 



Amn 2 {ui — ui) 



(1 + n 2 (ui - Mi)) 2 

Since v — z — AXf/n one can use the Jacobi elliptic functions addition theorem to obtain: 



cn iv 



z\m) 



cn(v 


m)cn(z 


m) + sn(v\m)sn(z\m)dn(v 


m)dn(2: 


m) 


1 — msn 2 {v 


m)sn 2 (z 


m) 



(18a) 
(18b) 
(18c) 
(18d) 
(18e) 
(18f) 

(19) 



This is a non-linear equation for a (i.e. for if> or u with respect to the type of the orbit). Equations 
ffT7al - [TBI except UM and FfSdl only give squares, so functions cn and sn are determined only 
up to the sign; dn is always positive. In order to determine those signs we plot in Figure 0] the 
functions sn and cn along orbits together with the interval where the orbit is defined. Figure 0] left 
shows that for type A orbits, only the function cn changes sign at the periastron and the function 
sn is always positive along the orbit. Thus, the sign of cn is positive if the orbit did not pass the 
periaston and negative otherwise. For orbits of type B (Figure |U middle) sn is always positive, 
while the sign of cn is determined by ()18aj) and (j!8d|) . thus no sign ambiguity arises. For orbits of 
type C the sign of cn is again unambiguous ( ([18af and (|18d|l ) while the function sn changes sign 
from negative before apastron to positive after apastron (Figure right). We conclude that the 
right hand side of equation ()19)) has two branches if it contains a sign ambiguity. The cases are: 
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Figure 5: Determining the angular momentum parameter for a type A orbit. Left: the orbit from (Pj to 
y f does not pass periastron. Right: the orbit from CPj to CPy- passes periastron. The designations right x 
and right 2 refer to the branches (|2U|) and (|2Tj) respectively. 

type A "right/' and "right 2 " and type C "right/' and "right 3 ", where 

cn(f \m)cn(z\m) + sn(t> \m)sn(z\m)dn(v \m)dn(z\m) 

nghtj = - tt— j r r-r-j 7 (20) 

1 — msn z {v\m)sn z {z\m) 

. —cn(v\m)cn(z\m) + sn(v\m)sn(z\m)dn(v\m)dn(z\m) . . 

nght 2 = - rr— ; 7 77— i 7 (21) 

1 — msm(t> |7?T.)sn^(^|mJ 

. cn(t>|m)cn(z|m) — sn(f |m)sn(^|m)dn(v|m)dn(2;|m) 

nght 3 = - r— : 7 7-7— j 7 • (22) 



In the above the sign of cn or sn is positive when calculated from a square root. 

Some examples of determining the value of a are shown in figures They show the left 
hand side and all the appropriate branches of the right hand side of (flU|) as a function of ip or 
ft for different cases of r iy rj and AA/. The solution of equation (TEH), which is the point where 
"right" crosses "left", is found numerically by using Brent method ((Press et al.lll988f ). 3 With the 
now known parameter a it is straightforward to calculate the value of to (equations llOJlrSjfKj) . 



3 Calculation of the time of flight 

The time of flight for a photon, is given by the integral (jChandrasekharlll992l b 



'/< = ±2Ma I -77, ; 7% v • ( 23 ) 



u 2 (l — u)\Ja 2 — u 2 (l — a 



The indefinite integral in the above formula can be expressed in analytic form using elliptic inte- 
grals. After considerable algebraic manipulations we obtain the analytic forms for the three cases 
of type A, B and C orbits as follows. 



3 For A A/ > 2ir there may be more than one solution to the equation : clearly in this case the orbit is wound 
about the black hole at r = 3M, therefore only the solution with a closest to ^= applies. 
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Figure 6: Connecting IP, and 7 f with a type B orbit and a type A winding orbit (winding number 
k = 1); no k = type ^4 orbit exists in this case. Left: the orbit from J", to CPy of type A. Right: the 
orbit from IPj to Tf of type -B . 




Figure 7: Determining the angular momentum parameter for a type B orbit. Left: a direct orbit from 
7i to 7f. Right: a winding orbit from 7i to 7 f (winding number k = 1). 
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Figure 8: Determining the angular momentum parameter for a type C orbit. Left: the orbit from CPj to 
Tf does not pass apastron. Right: the orbit from 3>j to passes apastron. The designations right} and 
right 3 refer to the branches ()2Uj) and (|22j) respectively. 




Figure 9: Connecting 5>j and Tf with a type B orbit and a type C winding orbit (winding number 
k = 1); no k = type C orbit exists in this case. Left: the orbit from T,, to of type C . Right: the 
orbit from to Tf of type B . 
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Type A: Introduce x such that fcf l9"fl) 

u = u 2 — (u 2 — u 3 ) cos 2 x 

to obtain: 



(24) 



t(x) 



2an 



lir. 



1+% + 



m/2 



2(m — ni)(r7 1 — 1) 



n(ni;x|m) + 



1 - 7i 3 



n(n 2 ; xM 



+ 7 T7 tt %H- 1 F(xm)-— — - 

(m-ni)(ni - 1) \ Vl ; V nj VAI ' 2(1 - m sin 2 x) 



m sm x 



(25) 



Here E, F and II are the incomplete elliptic integrals of the first, second and third kind respectively 
(jWolframl [l996) : u\, u 2 , U3, m and n are those from (|9*al - l9*e )l . and ni, n 2 are: 



1 u 2 

Tli = 1 

U3 
U 2 - M 3 



n 2 



I-U3 



(26) 
(27) 



Type B and C: Introduce x such that (cf 112 



u = U\ + 



1 1 — cos x 
n 2 1 + cos x 



to obtain: 



*(X) = 2a 



where 



ot 2 {n\ — 1) / ?7 1 sin x{ki cos x + 1) \/l — m sin 2 



nf + 777 — 277777-! 



fci(l — 77-1 sin x) 

1 / A; 
+ a 3 (l - 77 2 )(1 + — ) n(n 2 ; x\m) + - 2 In |x 2 | 
k 2 V 2y/n 2 - m 

1 / fci 
+ -ni)(l + — ) n(ni;x|»") + 7 In^il 

fel \ 2-^771 — 777 



- E(x\m) 



« 2 



+ F( X \m) -1(1 + (1 + fci)^(m - 1)) - t 1 " 



77 



Oi 

fcl 

m 

Xi 

:r 2 



n 2 u\ + 1 

1 — 7ii?7 2 
1 + 7Ji77 2 

kf-l 



a 2 
k 2 



77 



7? 



n 2 



{n 2 u\ + l) 2 

1 + (1 - Ui)77 2 
1 - (1 -7Ji)77 2 



«3 



77 2 (1 — Ui) — 1 



^'1 



\/t7i — 777 sin x + a/1 — 777 sin 2 



X 



y/ni — ?77 sin x — y 1 — ui sin x 



a/ ?72 — 777 sin x + y 1 — ?77 sin x 



y/n 2 — ?77 sin x — y 1 — "H7 sin x 



(28) 



(29) 



(30) 
(31) 
(32) 

(33) 
(34) 
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Figure 10: Graph r = r(t) in units of M for all three types of orbits. The type A orbit starts at infinity 
at t — ► — oo, passes the periastron at t = and continues to infinity as t — > oo. The type B orbit starts 
at the event horizon at t — > — oo, crosses the critical r = 3M at t = and continues to infinity as t — > oo. 
The type C orbit starts at the event horizon at t — > —oo, passes the apastron at t = M and again 
approaches the event horizon as t —* oo. In the above example a = — 10~ 3 for type A and C orbits 

and a = + 5 • 10~ 4 for type B orbit. 



For type B orbits the parameters Ui, m and n are those from (jl2al - [T2q) . while for type C orbits 
these parameters are from (|14al - IT3ej) . 

In figurelTTllwe show examples of r(t) for all three types of orbits. We compared the effectivness 
of the analytic method and the pure numerical method to calculate the times of flight. It was found 
that the algori thm based on the analytic method and implemented with the Carlson's algorithm 
( Carlson 1979j ) is always more accurate and 3 to 6 times fa ster than fou rth-order Runge-Kutta 
integration with adaptive stepsize control. 4 We also note (jKostic 2003) that it is possible to 
calculate the time of flight numerically by expanding the integrand (f2*3*j) into piecewise rapidly 
convergent series of analytically integrable functions. Such series easily give results accurate to 
10 _8 M for subcritical orbits and are some 6 times faster than the analytic method. The most 
effective method to calculate the time of flight would thus be a combination of the analytic method 
for the close to critical orbits and the series solution for the rest. 



4 Conclusions 

In this work we presented the complete solution of the ray-tracing problem in the Schwarzschild 
space-time. All the algorithms presented here have been tested for accuracy and for speed of 
execution and were found to be more accurate and considerably faster than commonly used direct 
integration methods. We hope that algorithms presented here will be used as a useful tool in 
solving more complex ray tracing problems that will elucidate the physics governing complicated 
transient phenomena in the vicinity of black holes. We would like to remind the community that 
a similar ray-tracing problem in the Kerr space-time still remains to be solved. 

4 We note that the algorithm ellpi of Numerical Recipes ijPress et al.lll98^l does not give the same resu lts as 
Mathematica for II integrals. We rewrote the function rj according to Carlson's original paper l|Carlsonlll979|) . and 
obtained identical results with Mathematica. 
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